/*
 * Basic component library: tapers
 *
 * Copyright (c) 1993-2010 Leo Spiekman, Xaveer Leijtens
 *
 * This file is part of libMask.
 *
 * libMask is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * libMask is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with libMask.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "mask.h"

#ifndef MIN
#define MIN(x,y)    ((x)<(y)?(x):(y))
//double MIN(double x, double y) {return ((x)<(y)?(x):(y));}
#endif
#ifndef MAX
#define MAX(x,y)    ((x)>(y)?(x):(y))
//double MAX(double x, double y) {return ((x)>(y)?(x):(y));}
#endif
#define SQR(x)      ((x)*(x))
//to avoid error due to SWIG.
//double SQR(double x) {return ((x)*(x));}

#define VERYSMALL   (1.0e-12)

double WMIN = 4e-6; /* Default contact lithography; use 1.0 for Canon */

static void taper_poly(double l, double w1, double w2)
{
  double x1,x2,x3,x4,y1,y2,y3,y4;

  movept(l,w1,4,7);
  x1=cp.x;
  y1=cp.y;
  movept(l,w1,7,1);
  x2=cp.x;
  y2=cp.y;
  movept(l,w1,1,6);
  movept(l,w2,6,3);
  x3=cp.x;
  y3=cp.y;
  movept(l,w2,3,9);
  x4=cp.x;
  y4=cp.y;
  movept(l,w2,9,6);

  mskpolygon4(clayer,x1,y1,x2,y2,x3,y3,x4,y4);
}

static void taper_angle_poly(double l, double a, double w1, double w2)
{
    double x1,x2,x3,x4,y1,y2,y3,y4,d;

    movept(l,w1,4,7);
    x1=cp.x;
    y1=cp.y;
    movept(l,w1,7,1);
    x2=cp.x;
    y2=cp.y;
    movept(l,w1,1,6);
    d=l*tan(a);
    movept(l,w2+2*d,6,3);
    x3=cp.x;
    y3=cp.y;
    movept(l,w2,3,9);
    x4=cp.x;
    y4=cp.y;
    movept(l,w2,9,6);

    mskpolygon4(clayer,x1,y1,x2,y2,x3,y3,x4,y4);
}

static void taper_rect(double l, double w1, double w2)
{
  double h,a,diag,ddiag,w;
  int i,j,n;

  h=(w2-w1)/2;
  a=atan2(h,l);
  ddiag=fabs(WMIN*h/l);
  diag=sqrt(l*l+h*h)+ddiag;
  if (w1>w2)
    ddiag=0.0;

  push_cp();
  movept(l,w1,4,7);
  move(0.0,0.0,a,NOFLIP);
  move(-ddiag,0.0,0.0,NOFLIP);
  rect(diag,WMIN,7,7);
  pop_cp();

  push_cp();
  movept(l,w1,4,1);
  move(0.0,0.0,-a,NOFLIP);
  move(-ddiag,0.0,0.0,NOFLIP);
  rect(diag,WMIN,1,1);
  pop_cp();

  n=(int)(fabs(h)/(0.5*WMIN*cos(a)))+1;
  for(i=0;i<n;i++) {
    j=(w2<w1)?n-1-i:i;
    w=MIN(w1,w2)-WMIN+2*j*fabs(h)/n;
    w=MAX(w,WMIN);
    sw(l/n,w);
  }
}

static void taper_angle_rect(double l, double a, double w1, double w2)
{
  double d,dl,dh,h1,h2,a1,a2,diag1,diag2,ddiag1,ddiag2,w;
  int i,j,n;
 

  d=l*tan(-a);
  h1=(a<0)?(w2-w1)/2+d:-(w1-w2)/2+d;
  h2=(a<0)?(w1-w2)/2+d:-(w2-w1)/2+d;
  a1=atan2(h1,l);
  a2=atan2(h2,l);
  
  ddiag1=fabs(WMIN*h1/l);
  ddiag2=fabs(WMIN*h2/l);
  diag1=sqrt(l*l+h1*h1)+ddiag1;
  diag2=sqrt(l*l+h2*h2)+ddiag2;
  if (w1<w2) {
    ddiag1=(a<0) ? ddiag1:0.0;
    ddiag2=(a<0) ? 0.0:ddiag2;
    dl=(a<0) ? cos(a1)*0.5*WMIN/(fabs(h1)/l):cos(-a2)*0.5*WMIN/(fabs(h2)/l);
  } else {
    ddiag1=(a<0) ? ddiag1:0.0;
    ddiag2=(a<0) ? 0.0:ddiag2;
    dl=(a<0) ? cos(a2)*0.5*WMIN/(fabs(h2)/l):cos(-a1)*0.5*WMIN/(fabs(h1)/l);
  }
  push_cp();
    movept(l,w1,4,7);
    move(0.0,0.0,a1,NOFLIP);
    move(-ddiag1,0.0,0.0,NOFLIP);
    rect(diag1,WMIN,7,7);
  pop_cp();
  push_cp();
    movept(l,w1,4,1);
    move(0.0,0.0,a2,NOFLIP);
    move(-ddiag2,0.0,0.0,NOFLIP);
    rect(diag2,WMIN,1,1);
  pop_cp();

 /* push_cp(); */
    n=(int)(l/dl)+1;
    move(0,d/l*dl/2,0,NOFLIP);
    for(i=0;i<n;i++) {
      j=(w2<w1)?n-1-i:i;
      w=MIN(w1,w2)-WMIN+2*j*fabs(w2-w1)/2/n - fabs(d)/l*dl;
      w=MAX(w,WMIN);
      sw(l/n,w);
      dh=d/n;
      move(0,dh,0,NOFLIP);
    }
  /* pop_cp(); */
}

void taper(double l, double w1, double w2)
{
  if (fabs(l) < VERYSMALL) return;
  if (mode&PATGEN) {
    WMIN = mode&CANON ? 1e-6 : 4e-6;
    taper_rect(l,w1,w2);
  } else
    taper_poly(l,w1,w2);
}

void taper_parabolic(double l, double w1, double w2)
{
  int READY=0;
  double y0, y1, y2, a;/* y=a*x*x */
  /* Maximum distance allowed between edges of lineair and parabolic taper.
   */
  const double dx = getlayer_accuracy(currentlayer());
  double w_tap1, w_tap2; /* start and end width lineair taper*/

  if (fabs(w1-w2)/(4*l) < VERYSMALL) {
    rect(l, (w1 + w2) / 2, 4, 6);
    return;
  }
  if (fabs(l) < VERYSMALL) return;
  if (w1>w2){
    skip(l);
    push_cp();
    rotate(M_PI);
    taper_parabolic(l, w2, w1);
    pop_cp();
  } else {
    a=4*l/(w2*w2-w1*w1);
    y0=a*(w1/2)*(w1/2);
    y1=y0;
    w_tap1=w1;
    while(!READY){
      w_tap2=w_tap1+4*dx+4*sqrt(dx*(w_tap1+dx));
      y2=a*(w_tap2/2)*(w_tap2/2);
      if (y2-y0>l){
	READY=1;
	y2=y0+l;
	w_tap2=w2;
      }
      taper(y2-y1, w_tap1, w_tap2);
      w_tap1=w_tap2;
      y1=y2;
    }
  }
}

void taper_parabolic_old(double l, double w1, double w2)
{
  int READY=0;
  int REVERSE=0;
  double temp;
  double y0, y1, y2, a;/* y=a*x*x */
  double dx=5e-8;  /* max. difference allowed between edges */
  /* lineair and parabolic taper. */
  double w_tap1, w_tap2; /* start and end width lineair taper*/

  if (fabs(w1-w2)/(4*l) < VERYSMALL) {
    return;
  }
  if (fabs(l) < VERYSMALL) return;
  if (w1>w2){
    REVERSE=1;
    skip(l);
    push_cp();
    rotate(M_PI);
    temp=w1;
    w1=w2;
    w2=temp;
  }
  a=4*l/(w2*w2-w1*w1);
  y0=a*(w1/2)*(w1/2);
  y1=y0;
  w_tap1=w1;
  while(!READY){
    w_tap2=w_tap1+4*dx+4*sqrt(dx*(w_tap1+dx));
    y2=a*(w_tap2/2)*(w_tap2/2);
    if (y2-y0>l){
      READY=1;
      y2=y0+l;
      w_tap2=w2;
    }
    taper(y2-y1, w_tap1, w_tap2);
    w_tap1=w_tap2;
    y1=y2;
  }
  if (REVERSE) pop_cp();
}

void taper_angle(double l, double a, double w1, double w2)
{
  if (fabs(l) < VERYSMALL) return;
  if (mode&PATGEN) {
    WMIN = mode&CANON ? 1e-6 : 4e-6;
    taper_angle_rect(l,a,w1,w2);
  } else
    taper_angle_poly(l,a,w1,w2);
}

/* Linear taper for two waveguides of length l, width w spaced
 * by gap g.
 */
void taper2(double l, double w1, double w2, double g1, double g2)
{

  double a;
  a = atan2((w2 + g2 - (w1 + g1)) / 2, l);

  push_cp();
  movept(l, w1 + g1, 4, 7);
  taper_angle(l, -a, w1, w2);
  pop_cp();
  push_cp();
  movept(l, w1 + g1, 4, 1);
  taper_angle(l, a, w1, w2);
  pop_cp();
  movept(l, w2 + g2, 4, 6);
}
